flist = dir("*description.mat");

temperature = [];
strain_data = [];
for j=1:length(flist)
    load(flist(j).name);

    strainnames = unique([sampleplate.description]);
    
    % Append temperature vector
    T = mean(sampleplate(1).temperature);    
    temperature = [temperature; T];
    

    G = [];
    E = [];
    N = [];
    for k=1:length(strainnames)
        sname = strainnames(k);
        I = [sampleplate.description]==sname;
        gmax = [sampleplate(I).max_growthrate];
        
        gmax_width = [sampleplate(I).max_growthrate_peakwidth];
        
        gmax_cut = gmax(gmax_width>8);
    
        gmax_mean = mean(gmax_cut);
        gmax_N = length(gmax_cut);
        gmax_err = std(gmax_cut);

        G = [G; gmax_mean];
        E = [E; gmax_err];
        N = [N; gmax_N];
    
    end

    strain_data = [strain_data G E N];

end


%% Plot growth vs. temperature across all strains and temperatures
close all;
figure;
colorlist = distinguishable_colors(length(strainnames));
for k=1:length(strainnames)

    gr = strain_data(k,1:3:end)';
    gr_err = strain_data(k,2:3:end)';
    errorbar(temperature, gr, gr_err, 'ko', 'markerfacecolor', colorlist(k,:));
    hold on;
%     plot(temperature, gr, 'color', colorlist(k,:))
end
box off;
set(gcf, 'Position', [0 0 400 300]);
set(gca,'fontsize', 20);
xlabel('Temperature (°C)');
ylabel('Growth rate (1/h)');
xlim([25 50]);
ylim([0.5 3.5]);
legend(strainnames);


%% Plot Arrhenius data and calculate activation energy
figure;
Ea_tot = [];
Ea_tot_err = [];

for k=1:length(strainnames)

    gr = strain_data(k,1:3:end)';
    gr_err = strain_data(k,2:3:end)';
    errorbar(1./(temperature + 273), log(gr), gr_err./gr, 'o', 'color', colorlist(k,:));
    tempcut = temperature<39;
    [f1,f2] = fit(1./(temperature(tempcut) + 273), log(gr(tempcut)), 'poly1');
    Ea = -f1.p1/298;
    ci = confint(f1);
    Ea_err = abs(Ea + ci(1,1)/298);


    % Compute error in fit coefficients
    alpha = 0.95; % 95% confidence interval
    ci = confint(f1,alpha); % Obtain intervals of confidence (95%)
    df = sum(tempcut); % Degrees of freedom in T-test
    t = tinv((1+alpha)/2, df); % Use student's T test distribution
    se = (ci(2,:)-ci(1,:)) ./ (2*t); % Standard error of T-test
    
    Ea_tot = [Ea_tot; Ea];
%     Ea_tot_err = [Ea_tot_err; Ea_err];
    Ea_tot_err = [Ea_tot_err; se(1)/298];

    hold on;
    plot(1./(temperature(tempcut) + 273), f1(1./(temperature(tempcut) + 273)), 'color', colorlist(k,:));
    
end
box off;
set(gcf, 'Position', [0 0 400 300]);
set(gca,'fontsize', 20);
xlabel('1/T (1/K)');
ylabel('Log(growth rate)');
xlim([25 50]);
ylim(log([0.5 3.5]))


%% Version 2 of Arrhenius plot
figure;
Ea_tot = [];
Ea_tot_err = [];

for k=1:length(strainnames)

    gr = strain_data(k,1:3:end)';
    gr_err = strain_data(k,2:3:end)';
    plot(1./(temperature + 273), log(gr), 'ko', 'markerfacecolor', colorlist(k,:));
    tempcut = temperature<39;
    [f1,f2] = fit(1./(temperature(tempcut) + 273), log(gr(tempcut)), 'poly1');
    Ea = -f1.p1/298;
    ci = confint(f1);
    Ea_err = abs(Ea + ci(1,1)/298);


    % Compute error in fit coefficients
    alpha = 0.95; % 95% confidence interval
    ci = confint(f1,alpha); % Obtain intervals of confidence (95%)
    df = sum(tempcut); % Degrees of freedom in T-test
    t = tinv((1+alpha)/2, df); % Use student's T test distribution
    se = (ci(2,:)-ci(1,:)) ./ (2*t); % Standard error of T-test
    
    Ea_tot = [Ea_tot; Ea];
%     Ea_tot_err = [Ea_tot_err; Ea_err];
    Ea_tot_err = [Ea_tot_err; se(1)/298];

    hold on;
    plot(1./(temperature(tempcut) + 273), f1(1./(temperature(tempcut) + 273)), 'color', colorlist(k,:), 'linewidth', 1);
    
end
box off;
set(gcf, 'Position', [0 0 400 300]);
set(gca,'fontsize', 20);
xlabel('1/T (1/K)');
ylabel('Log(growth rate)');
%xlim([3.17E-3 3.35E-3]);
xlim([3.21E-3 3.35E-3]);
ylim(log([1 3.5]))

%% Plot bars with errors of activation energies (in units of kcal/mol)
x = 1:length(strainnames);
figure;
% bar(x,Ea_tot);                
hold on
for k = 1:length(strainnames)
    er = errorbar(x(k),Ea_tot(k)/1.688,0,Ea_tot_err(k)/1.688, 'markeredgecolor', colorlist(k,:)); 
    bar(x(k),Ea_tot(k)/1.688, 'facecolor', colorlist(k,:));
    hold on;
%     er.MarkerFaceColor = colorlist(k,:);
    er.Color = [0 0 0];                            
%     er.LineStyle = 'none';  
end

box off;
set(gcf, 'Position', [0 0 400 150]);
set(gca,'fontsize', 20);
ylabel('Activation energy');
ylim([0 20])

